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Abstract: We analyze the modulation instability induced by periodic 
variations of group velocity dispersion and nonlinearity in optical fibers, 
which may be interpreted as an analogue of the well-known parametric 
resonance in mechanics. We derive accurate analytical estimates of resonant 
detuning, maximum gain and instability margins, significantly improving 
on previous literature on the subject. We also design a periodically tapered 
photonic crystal fiber, in order to achieve narrow instability sidebands at a 
detuning of 35 THz, above the Raman maximum gain peak of fused silica. 
The wide tunability of the resonant peaks by variations of the tapering 
period and depth will allow to implement sources of correlated photon pairs 
which are far-detuned from the input pump wavelength, with important 
applications in quantum optics. 
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1. Introduction and motivations 

Parametric resonance (PR) is a well-known instability phenomenon which occurs in systems 
the parameters of which are varied periodically during evolution, see the classical treatments in 
Refs. (T3I2J. For example, a harmonic oscillator the frequency of which is forced to vary in time, 
will become unstable if its internal parameters and the amplitude of the frequency variation 
happen to be inside special regions, known as resonance tongues. The study of the properties 
of resonance tongues has a long history and relies on a variety of geometrical approaches E|4]. 

It is natural that such a general phenomenon was associated to the equally important instabil- 
ity process that is ubiquitous in infinite dimensional dynamical systems: modulation instability 
(MI), also known as Benjamin-Feir instability EIS). MI is known to exist in different branches 
of physics such as fluid-dynamics (7|, plasma physics @[9][10), Bose-Einstein condensates 
ifTTl and solid-state physics |[T2l . In nonlinear optics 1131 . it manifests itself as an exponential 
growth of two spectrally symmetric sidebands on top of a plane wave initial condition, by virtue 
of the interplay between the cubic Kerr nonlinearity and the group velocity dispersion (GVD). 
In optical fibers it leads to the breakup of a plane wave into a train of normal modes of the 
system, i.e. solitons lfl4l[l5l . 

The link between PR and MI has been established in 1993, see lfl6l . in relation to the periodic 
re-amplification of signals in long-haul telecommunication optical fiber cables. This was based 
on a nonlinear Schrodinger equation (NLS) where the coefficient of the nonlinear term is varied 
along the propagation direction. Importantly, this peculiar type of MI occurs in both normal and 
anomalous GVD. This prediction was later partially verified in experiments, see lfT7l[T8l . 

Moreover, in long-haul fibers, dispersion management is a commonly used technique which 
introduces periodic modulation of fiber characteristics. The possibility of instability phenomena 
disrupting adjacent communication channels has been thoroughly analyzed, see e.g. IU9ll20ll2Tl 
l22l . Specifically, in lfl9l the partial suppression of the conventional MI in anomalous GVD due 
to a large swing dispersion management is discussed, while in l20l the degenerate case of zero 
average dispersion was studied. The combination of both loss and dispersion compensation is 
studied in 12111221 . The main interest in those works was on step-like variations of the GVD 
coefficient. 

At the same time the effects of smooth periodic or random variations of fiber parameters 
were studied in [23. 24. 25]. Also some work has been done on the effect of the perturbation of 
fiber parameters on soli ton propagation [26. 1271 . 

It turns out that the variation of dispersion and nonlinearity can enhance or suppress the 
PR, while higher order nonlinear effects such as self-steepening proved less important. Quite 
surprisingly, experiments on micro-structured fibers have been reported only recently for the 
first time, see Refs. SUED, where a photonic-crystal fiber (PCF, ED) of varying diameter 
is used. In that experiment, the dispersion is periodically switched from normal to anomalous, 
but this feature is not required to achieve PR, while the effect of Raman scattering plays an 
important role in the relative magnitude of the PR peaks. 

The conventional explanation is in term of a grating-assisted phase matching process 1161 
l2Tll22ll28l . with the exception of Ref. Il24l . which relies on the theory of Mathieu equation — a 
paradigmatic example of the application of Floquet theory, see Refs. (2] O — and analyzes the 
case of nonlinearity variations only. 

In this work we present an improved analytical treatment based on the so-called Poincare- 
Lindstedt perturbation method and averaging method, which provide excellent predictions on 
resonant frequency, gain and bandwidth of the unstable peaks appearing in the output spectra. 
The new feature here is the tunability of the instability bands, specifically the lowest order one, 
which possesses a gain larger than the other PR peaks. Finally we provide physical estimates 
for a PCF-based system which leads to instability bands with frequency detunings larger than 



13 THz, the detuning of the Raman gain peak of fused silica. This system can be used for 
the generation of entangled photon pairs in a narrow bandwidth, far detuned from the Raman 
gain peak, which would reduce the impact of Raman-induced decorrelations, with important 
applications in quantum optics, especially in quantum computation and cryptography II3T1I321 . 

2. Model equation and Floquet theory: analytical and numerical analysis 

2.1. NLS with varying parameters and derivation of Hill equation 

Let us consider the dimensionless NLS equation for a slowly-varying electric field envelope 
A(z,t) of a guided mode at carrier frequency ©o, with both GVD and nonlinearity varying 
along the propagation direction, namely 

id-A - l -s{z)d?A + n(z)\A\ 2 A = 0. (1) 

s and n are normalized coefficients, s(z) = & (z) /P® an d n ( z ) = Y( z ) I T°> wnere Pi (z) an d /(z) 
are the GVD and nonlinear coefficients, respectively, and the superscript denotes their mean 
values, n and s are assumed to be periodic functions of z. Finally z = Z/Z„; is the dimensionless 
distance in units of the nonlinear length Z„/ = (y /^) -1 , and t = (T — v Q g Z)/T s is the dimension- 
less retarded time in units of T s = yZ^i\P^\, and v° is the mean group velocity. P t is the total 

input power injected in the mode, and A is the dimensionless modal intensity scaled by y/P~ t . 

We look for for a steady state solution of (HJ in the form A = \/Pq exp (i(j> (z)) : it can be ver- 
ified that 0(z) = Pq Jz oo n(z / )dz' . We then perturb this steady state by adding a small complex 
time dependent contribution a(z,t), i.e. A(z,t) = (\/^) + £a(z,t)) exp(/0(z)), with £< 1. In- 
serting this ansatz in Eq. ([]]) and taking only the terms which are first order in e, one finds that 
a obeys the following equation: 

1 T 

id z a - -s(z)d?a + n(z)P (a + a*) = 0. (2) 
Finally we substitute in (O the ansatz 

a(z,t)=a A (z)e- iSt + a s (z)e iSt 
to obtain the following Schrodinger equation 

=H s (z)\y), H s (z) = (-j(z)y -«(z)/b) 6 z -in(z)Po6 y (3) 

where the dot denotes z-derivative, H s (z) is a non-Hermitian Hamiltonian (which thus allow 
unstable modes to grow), \ \jf) = {aA,a* s ) J ', and d,- are the Pauli matrices. 

In the remaining part of this section we will assume that dispersion and nonlinearity exhibit 
the simplest possible periodic behavior 

s(z) = so + s(z) = sq + hs\ cos az, n(z) = «o + «(z) = »o + hn\cos az, (4) 

where generally so = ± 1 for normal (anomalous) dispersion and no = 1 ; oc is the normalized 
spatial angular frequency for the parameter oscillation. The forcing amplitude is controlled by 
parameter h, which must be small to guarantee the validity of our perturbative expansions. 
However, we find below that our estimates are reliable even for h ~ 0.5. 

The conventional way in which PR is approached in fiber optics is to split the dispersion co- 
efficient into two parts: a constant average term and an oscillating term. By means of a change 



of variables this latter is included in a single varying nonlinear coefficient, see Refs. 1191 1211 . 
Then all quantities are expanded in Fourier series and it is assumed that only one Fourier coef- 
ficient resonates at each PR order. This leads to the nonlinear phase-matching condition for the 
rn-th peak resonant frequency 

s 8l+2n Po = ma (5) 

where m is the PR order and S m the corresponding resonance detuning. This condition is valid 
for m>0 (m < 0) if GVD is normal (anomalous). Below, in our derivation, we assume instead 
that m is a positive integer. Eq. (|5]l can also be understood as a grating assisted phase matching 
condition, the grating being the periodic modulation of the fiber parameters. We performed 
detailed numerical simulations and verified that the result of (O is quite coarse. The purpose 
of the present section is to prove that a more accurate estimate can be made, in a more general 
case than the variation of nonlinear coefficient only, which was already discussed in 1241 . We 
define the parameter space as the plane (S,h) and we investigate the different features of PR 
instability (resonant frequency, bandwidth and gain) on this plane, for a fixed value of \n\ \/\s\ \. 

In order to compare to the theory of PR in a classical mechanical oscillators, it is instructive to 
derive a 2nd order ODE from the system of Eqs. (O. Let us first transform it in phase-quadrature 
variables, by applying the rotation 



We define the new state |0) =R\\j/) which satisfies a system analogous to that in Eq. (O, with 
Hamiltonian matrix 

Cl (z) N 



with c\ (z) = — js(z)8 2 and C2(z) — c\ (z) — 2«(z)Po. We can thus more easily derive a second 
order ODE for both 01,2: 

C\ 2 ■ 

01.2 -01,2 + C1C 2 01,2 =0, (7) 

Cl,2 

which still contains a first derivative of 0i,2- This can be eliminated by the substitution 

01.2 = exp (\ ( C ^ 1,2 dz /N ) 01.2 = VCL201.2- 

\2Jq c(z') 1,2 / 



We finally obtain the following ODE 

01.2+<C 1C 2 + - 

2ci 2 4 



• n 2" 

Cl,2 



c l,2 



01,2 = (8) 



Notice that in contrast to a usual dissipative term, which corresponds to a threshold in the 
magnitude of the perturbation needed to achieve PR instability (see Ref. |f2]), our variable 
transformation is simply an oscillating factor, which does not affect the instability. 

The left-hand side of (O contains, even in our simple case of Eq. ©, several harmonic terms. 
This means that Eq. (O is a Hill equation, which generalizes the Mathieu equation found in [24 1, 
where c\ = 0. It is well known that Mathieu equation is an exceptionally simple case, HHU, 
with a predictable structure of stability and instability regions in the parameter space. Moreover 
since we consider a Hill equation, we expect the instability regions to be irregular and appear in 
the form of instability islands separated by stable parts. This is consistent with the statements 
made in Ref. |[T9l . i.e. that a large switching of dispersion may suppress instability. 



Despite Eq. ([H) highlights all these important properties of the system, it is difficult to handle 
analytically. We thus turn back to (0 and derive our analytical estimates from it, in order to 
report the simplest possible treatment. 

In the following part we present the main results of our calculations, which rely on the com- 
monly used methods of averaging and on the Poincare-Lindstedt perturbation method, see 

2.2. Estimate of Resonant Frequency and Parametric Gain 
We now explore the properties of the system 



i\j))=H pq (zM 



(9) 



with 7/pq defined in (O. In the limit of vanishing perturbation we have a simple harmonic 
oscillator written in complex variables 0i 2. It is widely known that, in this limit, parametric 
resonance occurs if the natural frequency of the oscillator is a multiple of half the forcing 
frequency and this is the basic point behind any perturbative approach, see J2J. Thus, by setting 

c i,2(z) = c ? 2 + ^12(2), the natural frequency of the unforced oscillator is simply (Oq = <J CjCj. 



The resonance condition becomes Oq = ma/2, with m = 1, 2, 3, . 
of the NLS parameters corresponds to a detuning 



which expressed in terms 



1 



2n P Q 



I { ma V 
+ \2noFa) 



(10) 



Note that, if we assume j^j^ ^> 1, we obtain the quasi-phase matching condition given in 
Eq. © (taking care of the different convention on m, which in our relation, Eq. [10] is only 
positive, while in Eq. (0 can also be negative), which is thus a coarse approximation if a « 
n P . 

In order to obtain the peak gain at resonance we use the method of averaging, which consists 
in posing 



0i = A(z) cos (Oqz + B(z) sin coqz, 02 = — — g- [A(z) sin coqz - B(z) cos coqz] 



(ID 



substituting in (O and averaging over a period of the forcing term T z = 2n/a, we obtain 
A\ a(p 2 -pr)coi f-l+cos^) sin UA 



BJ 2%{a 2 -Acol) 



sin 



1- COB 



(12) 



where pi = hs\/so and P2 = h(8 2 s\ +4«iPo) / {S 2 sq +4-hqPq). It is thus apparent that the m = 1 
resonance occurs at 2 too = a and in this case the matrix of the system ( fT2l ) has purely real 
eigenvalues of opposite sign. The positive one corresponds to the peak gain of the first PR band 
and turns out to be the maximum achievable value. This reads as 



2a 



h\s m 



SittO 



(13) 



where Si is calculated according to Eq. ( 110| >. 

We notice that there may exist a set of parameter values which suppresses or even forbids 
the occurrence of the instability, specifically s\/sq = «i/«o- This implies that, e.g., in normal 



GVD (so = 1) if both nonlinearity and dispersion undergo parallel increase and decrease, the 
instability is suppressed. Instead the same condition maximizes the gain under anomalous GVD 
fa) = -l). 

The higher order resonances are more difficult to characterize by this method, but their gain 
is generally of order h 2 , since in this case the first-order contribution vanishes, see Eq. ( 112b . 

2.3. Estimate of resonance bandwidth 

Floquet theory predicts that our system © has quasi-periodic solutions (composed by a pe- 
riodic function with period T- and a phase-factor, a complex function of unit absolute value). 
Moreover, stability margins correspond to a pair of periodic or anti-periodic solutions, defined 
by \<p(z + T z )) = ±|0(z)), which possess periods equal to T z and 2T Z respectively. Bearing this 
in mind, we apply the Poincare-Lindstedt method, by making the following perturbation ex- 
pansion in powers of h: 

- ^S 2 = d +hdi+h 2 d 2 , |0(z)) = |0o(z)) +/#i(z)) +h 2 Mz)) + ■■■, 10,) = Oi;, <ki) T ■ 



(14) 



We choose do such that it corresponds to the first resonant frequency, i.e. do = ~8f/2 
At zeroth order we obtain 

0io + C0 o 2 0io = O (15) 

The first resonant band occurs at Oo = a/2, thus Eq. (fTBT l has solutions with period 2T Z . For 
a fixed value of h, the stability margins correspond in general to different frequencies, which 
implies they must correspond to linearly independent eigenfunctions. The higher order correc- 
tions, which provide the stability margins, are obtained by solving for the successive terms of 
the perturbation series (|0i (z)) , | 02 (z) ),•••) by making use alternatively each of the independent 
solutions of Eq. ( 115l >. 

We thus first pose 0io = cos C0qZ, and at first order in h we obtain: 



«o0n 



- dx s Q (cUcl) + (cl-c\) d -^ + c\n lP » 



COS COQZ+ 

cos3c0oz. 



(3c° + c ( D d -^+c%P Q 

We then impose that the secular term vanishes and solve for d\, 

d i = ^ — 7T~ —d Po{s Q ni-nosi). (16) 

2s (doSQ-noPo) 

If we set 0io = sin (Oqz we obtain the same value with opposite sign, so that the instability 
margins of the first band satisfy 

8 2 8 2 

Y = iTh dl (17) 

The four solutions of Eq. (fTTI i are denoted by ±5j% and the bandwidth (i.e. the range of un- 
stable detuning values) is given by F\ = 8^ — <5f . Higher order corrections and bandwidth of 
overtone resonances can be found by solving for the higher order terms |0i (z)), |02(z)), ■ ■ ■ of 
the perturbative expansion, but the calculation is too lengthy to report here. 

2.4. Comparison of analytical and numerical results of the Hill equation 

The Hill equation ((HJ) can be solved numerically by means of an ODE solver, following the 
prescriptions of the Floquet theory. The monodromy matrix is calculated by setting two linearly 
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Fig. 1. Resonant tongues in the (8,h) plane for the parametric resonance instability of a 
NLS with varying dispersion and nonlinearity, for (a) the first and (b) the second resonant 
bands. GVD is normal (so = +1) and a = 10, n\ = —s\ = 1 are chosen in order to obtain 
the maximum gain and bandwidthm, see Eqs. dl3l > and (16). The colormap provides the 
values of gain computed by means of an ODE solver. The solid green lines are the analytic 
predictions of Eq. (XT) and the blue dotted lines with markers are the band edges computed 
numerically by employing the Hill determinant method (see in the text). The insets show 
the maximum gain vs. h (blue dotted line with markers) and in (a) also the curve obtained 
from Eq. d 1 3 1 > (green solid). Notice that since we reported in the text only the perturbation 
results of the first PR tongue, in (b) we show only numerical results. 

independent initial conditions and evaluating the solution at z = T z , see The eigenvalues of 
the monodromy matrix provide the instability gain (Floquet exponents). 

Since PR bands are typically narrow, this requires a fine grid in the (8,h) parameter space. In 
order to speed-up the calculation we (i) compute the exact instability margins, then (ii) calculate 
the gain values of each instability tongue. 

The calculation of the exact stability margins can be carried out by the standard Hill de- 
terminant method or harmonic balance based on the truncated Fourier expansion of the vari- 
ables (0i 2) and forcing terms, see 1341 . In contrast with more conventional spectral problems, 
where the eigenvalue appears explicitly in the equation, we need to find the values of detuning 
5 at a fixed h, while the coefficients depend on 8 in a nontrivial way. This implicit depen- 
dence is solved by means of a root finding algorithm based on the minimization of the least 
singular value, see 1351 . The second step involve the numerical evaluation of the Floquet ex- 
ponents in the close proximity of each band. We compare the results of the numerical calcula- 
tion of resonant tongues and analytical estimates for the case of maximal gain and bandwidth, 
i.e. «i/«o = si/so, see d 1 3b and ( 1161 ). 

First we show in Figs. Q] and [2] the instability regions for both normal and anomalous GVD 
at a fixed frequency of parameter variation a = 10 for the first two PR bands, m = 1 (a) and 
m = 2 (b). Each resonance tongue stems exactly from the detuning predicted by Eq. (TTOb and 
the maximum gain (shown in the insets) only slightly drifts away from that value. Our ana- 
lytical estimates refer to the first resonance band and are shown in Figs. 02a) and [2 a) to be 
quite accurate. Instead in Figs.QJb) and|2jb) only numerical results are provided. We observe 
that the second order PR exhibits a threshold value for h below which the instability gain is 
virtually zero. This occurs also for higher-order PR and both in normal and anomalous GVD. 
As explained above this is not due to the first derivative in Eq. (O which is not a damping term, 
but can be qualitatively ascribed to the the fact that in the Hill equation the forcing function 
contains several overtones of a and they could in special cases suppress completely a subset of 
resonances; see 031 and references therein. 




Fig. 3. Characterization of the first three PR peaks (m = 1,2,3) as a function of a. Nor- 
mal GVD (sq = +1), s\ = — 1, ti\ = 1 and h = 0.5. We plot in (a) the resonant frequency 
calculated in Eq. d 1 Ob and obtained from the ODE solution as the point which maximizes 
the instability gain, which is shown, in logarithmic scale, in (b), which includes also the 
solution of Eq. dl3l >. for m = 1. Finally in (c) the instability bandwidth calculated numer- 
ically by means of the Hill determinant method and analytically for m = 1 , according to 
Eq. (16). In every panel the same convention is used, i.e. solid lines represent analytical re- 
sults, while dotted lines with markers are obtained by numerical calculations: specifically 
the blue dotted line with crosses is associated to m = 1, the green dotted line with stars 
to m = 2 and the red dotted line with circles to m = 3. The solid lines use the same color 
convention. 



We provide also the curves of resonant frequency, gain and bandwidth as a function of a 
for fixed amplitude of the variation of nonlinearity and dispersion, \s\\ = «i = 1 and h = 0.5, 
Figs.[3]and|4] The forcing terms s\ and n\ are of equal amplitude as above, and the perturbation 
is quite large. Nevertheless our estimates for the first band prove quite reliable; moreover the 
resonant frequency hardly differs from its h — > value at any resonant peak. 

We notice that in the normal GVD regime (Fig. [3]) each resonant frequency converges to zero 
as a —> while the gain g —> h'" from below as a —> °°. Anomalous dispersion (Fig. [4j causes 
gain to approach the same limit but from above, while the resonant frequency is limited by the 
conventional MI band which cuts off at 8 = 2. 
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Fig. 4. Same of Fig.[5]in the case of anomalous GVD (sq =— 1) and 5i = +1. 



2.5. Including higher order effects 

To conclude this section we briefly discuss how to compute the resonant frequency in a gener- 
alized NLS model including higher order terms. As above the linearized equation can be cast 
as 

i\^)=H gois (z)\(l)), H gnls = H°+H(z). (18) 

where we split as before the Hamiltonian in two parts, a constant and an oscillating one. We also 
assume that any common diagonal term of H gn \ s is eliminated by a suitable change of variable, 
which is always possible. At that point we have a traceless 2x2 matrix, thus if we define ±©o 
to be the eigenvalues of H°, and assume they are both real (for a certain choice of parameters 

and detuning), the resonance condition becomes (Oq = \J c\c^ = a/2 and from the expression 
of C0q we derive the corresponding detuning. 

We consider a generalized model which includes higher-order dispersion (HOD) terms up to 
the fourth order. It is well known that only even HOD terms contribute to MI ll36l . We thus just 
need to modify ci 2 in Eq. $6^ by making the following substitution 

- l -s( z )S 2 ^-\s{z)S 2 -^Z8\ 

where is the normalized fourth order dispersion (FOD) and is assumed to be constant. The 
third order dispersion (TOD) appears only as a common diagonal term and can be removed. 
The resonance condition is thus a quartic equation in 8 2 and if < there may exist more 
than just one pair of sidebands which undergo parametric resonance. This is analogous to the 
conventional MI in the presence of HOD, see e.g. (37J. 

3. Comparison with numerical solution of NLS 

In order to assess the correctness of our analysis we solve Eq. ((TJ by means of the split-step 
Fourier method. We use the same parameters as above: sq = I, a = 10, — S\ = n\ = 1 and 
h = 0.5. We use a large a in order to achieve widely separated PR resonances and large gain in 
the normal GVD regime, see Fig. [3] Finally we operate under normal GVD, because anomalous 
GVD gives rise to the conventional MI bands which have a larger gain (four times larger than 
the PR gain in the case of /; = 0.5). 



In Fig. [3 a) we show the output spectrum after a propagation distance z = 38. One can clearly 
distinguish the first three peaks (m =1,2, 3), plus two peaks corresponding to the four-wave 
mixing of the carrier and the m = 1 PR band. In table [T] we report the values of peak detuning 
and bandwidth extracted from the spectrum of Fig. [2a), the numerical results of the linearized 
problem discussed in the previous section and predictions obtained by phase-matching consid- 
erations, Eq. ©. The peak positions are in very good agreement with the results of the previous 
section. The approximate phase-matching always underestimates the value of the actual detun- 
ing. 

In Fig. [gib) we show how the three PR peaks evolve during propagation: the instability leads 
to amplification at the rate predicted by the linearized problem, i.e. Eq. ( fT3l ) (see the dashed 
lines), apart from saturation occurring towards the end of the evolution. 

Next we discuss the small scale oscillations reported in 1281 [29), i.e. an amplification- 
deamplification cycle undergone by the resonant peaks. Since Fig. [2b ) does not allow us to 
visualize them, it is useful to introduce a new set of parameters which gives larger gain, so that 
fewer periods are sufficient. Thus we choose h = 0.9 and normal GVD conditions. Moreover, 
in order to accurately visualize the oscillations, we set a larger period, i.e. a = 5. In Fig. [3c) 
we show the evolution of the first peak, together with the solution of the averaged problem, 
Eq. (II lb . with coq = a/2 (a phase shift is applied so that the two curves overlap). The numer- 
ical and approximate solutions of the linear problem agree quite well, thus providing a good 
explanation of the process without explicitely resorting to a three wave dynamical system (as 
in ||29l ). The amplification-deamplification depends in practice exclusively on a, i.e. it stems 
from the forced oscillations impressed on resonance by the variation of parameters, superim- 
posed to the unstable growth. This is completely analogous to a pendulum the pivot of which 
is displaced periodically: a small displacement from the position in which the bob points down 
initiates oscillations (at the natural frequency) which are further amplified at each cycle. 

For the sake of completeness we compare, in table [2] the position of resonant peaks with 
the predictions of the grating phase matching and the Hill equation. Our estimate Eq. < l 1 b is 
the closest to numerical NLS result despite it is of zeroth order in h, while the simple phase- 
matching argument, Eq. (O, is evidently less and less accurate as a — >• 0. 



m 


8m 


r m 




NLS Hill ([T0]l Grating © 


NLS Hill (num.) 


1 


2.8463 2.8632 2.8284 


0.34558 0.27481 


2 


4.2537 4.2544 4.2426 


0.087965 0.081239 


3 


5.303 5.2978 5.2915 


0.050265 0.035274 



Table 1. Values of resonant frequencies and bandwidth for a = 10. Comparison between 
data extracted from NLS evolution and the results of linearized model (solved both as 
phase matching condition and Hill equation). The resonant frequency is estimated at any 
order analytically, while numerical values are reported for bandwidth. 



m 


8m 




NLS Hill dl(J Grating <|5j 


1 


1.8912 1.8399 1.7321 


2 


2.9719 2.8632 2.8284 


3 


3.6191 3.6239 3.6056 



Table 2. Same as in Tab.[T] but for a = 5. 




Fig. 5. Numerical evolution of the NLS equation: (a) comparison of input (blue dashed 
line) and output (black solid) spectra. GVD is normal and the total propagation length is 
z = 38, a = 10, h = 0.5 and n\ = —s\ = 1. The m = 1, 2, 3 PR peaks as well as the first 
two mixing products of m = 1 are highlighted by arrows, (b) Amplification of the first 
three peaks: extracted from the spectrum (solid line) and predicted by the linearized anal- 
ysis, i.e. exponential growth with gain g m (dashed line), (c) The detail of the amplification 
process on a shorter scale: we use a different set of parameters, a = 5 and h = 0.9 (with 
n\ = — S\ = 1 as above) in order to have higher gain and a larger period. Blue solid line 
represents the evolution of the 1st PR peak spectral component, the dashed green line the 
solution of the averaged equations, Eq. i ll It . The amplification-deamplification process is 
apparent and agrees with the prediction of the averaged linear equation. 




Fig. 6. Main parameters of a PCF as a function of pitch A at wavelength Ao = 1064 nm. 
The fiber is made of pure silica and the filling fraction d/A = 0.4. The red dashed lines in 
(a,b) show the range spanned by the parameters when tapering the fiber. 



4. The design of a periodically tapered PCF 

In this last section we propose a feasible system which supports PR instability peaks exhibit- 
ing relatively large frequency detunings. This is important, e.g., in quantum optics where the 
implementation of new sources of entangled photons with a reduced Raman decorrelation is 
of great practical interest, and this can provide an alternative approach to what proposed in 
Refs. Il3ni32l. 

In our calculations we use a periodically tapered PCF whose index contrast, small core and 
design flexibility allow to obtain a large nonlinearity, together with regions of small GVD. A 
similar system was already used in |28 , 29], but here we predict the possibility of observing far 
detuned, tunable instability peaks, by employing short tapering periods. 

At first we explore the design space (pitch and filling fraction) in order to operate in a region 
of small normal GVD near Ao = 1064 nm, and such that the zero-dispersion point (ZDP) can be 
approached by slightly varying the PCF geometry. The modal analysis is performed by means 
of COMSOL Multiphysics |38l 

In Fig. [6]we show the properties of a PCF made of pure silica with a triangular lattice of air 
holes. We assume the air filling fraction c//A = 0.4 {d is the air hole diameter, A is the pitch) 
to be constant so that by varying A, d is adjusted accordingly. We are interested mainly in two 
quantities: the GVD (fa), see Fig.|6ja), and the nonlinear coefficient (7), Fig.[6{b), calculated 
according to the full vectorial model of the effective area, see ||39"1 . Fig.|6]also shows the next 
two terms of the Taylor expansion of the GVD, j33 [Fig. [6jc)] and /it [Fig.|6jd)], as functions 
of the pitch A. 

We propose to operate between A m j n = 3.25 jJ.m and A max = 3.6 jim. In this range, 
p 2 6 [0.4,2.7] ps 2 /km and ye [7.8,9.4] W^'/km. Thus the GVD undergoes, when compared to 
its average value, very large oscillations, equivalent to h « 0.85 in Eq. ©. Note that we do not 
cross into the negative GVD region in order to inhibit any spurious occurrence of the classical 
MI. Moreover the dependence of GVD on A is, in good approximation, linear. Thus a cosine- 
shaped tapering of the PCF leads to a cosine variation of GVD. The value of 7 is only slightly 
modified by the tapering and can be also approximated as a cosine. This allows to straightfor- 
wardly apply the theory developed in the previous sections. We thus use the following simple 



formulas for the variations of the parameters: 

A = A + A 1 cos az, Pi = P® + Pi cos y = y° + y 1 cos az, 

where all the parameters are reported in table [3] We use as above the superscripts and 1 to 
denote the average value and first Fourier coefficient, and we use a for the dimensional spatial 
frequency of tapering (a = OC/Znl), As far as HOD terms are concerned, they are in a good 
approximation constant and are also reported in table [3] Finally in our simulations we include 
the self-steepening term and stimulated Raman scattering response of silica, see [36], in order 
to obtain a very realistic simulation. 



d/A 


0.4 


P, 


25 W 


A 


3.425 jum 


A 1 


0.175 ftm 


P\ 


1.3695 ps 2 /km 


PI 


1.1275 ps 2 /km 


f 


8.7302 /(W km) 


f 


-0.8 140 /(W km) 


ft 


6.43 x 10" 2 ps 3 /km 


Pa 


-0.97 x 10" 4 ps 4 /km 


a 


47.64 m" 1 


^taper 


225 m 


itot 


250 m 


number of periods 


1708 



Table 3. PCF parameters used in the generalized NLS model. 

We set our design target: the period of variation along the direction of propagation (Tz = 
2n/a) has to be such that the m = 1 PR band is located at Af = ±35 THz. This quite large 
detuning allows to clearly distinguish PR from Raman gain, in contrast with [28 1 where PR 
peaks occur in the Raman gain spectral region and thus are further amplified by it (in the cited 
paper the spectrum exhibits a frequency asymmetry which is a clue to the enhancement of 
sidebands provided by Raman amplification). 

We verified that self-steepening and Raman effects are of little importance in our case, never- 
theless at large detuning we cannot neglect the effect of HOD. As discussed above, in sec. 12.51 
only even order terms contribute to the instability, and amount to a simple modification of 
ci 2 in Eq. ©. In order to complete our design, we thus set P, = 25 W, correct the PR con- 
dition including f$4 and finally obtain a = 47.7 m -1 , corresponding to a period of tapering of 
Tz = 13.2 cm. Neglecting ($4, the estimated period would be less than 10 cm. The predicted 
gain at Af = ±35 THz is approximately g w 100 km -1 . We simulate a 250 m long fiber with a 
periodically tapered central part of about 225 m, which corresponds to 1708 periods. The power 
level and fiber tapering periods are only slightly more demanding than those used in 11281 . The 
use of a highly nonlinear fiber (for example by using materials other than fused silica) can scale 
down power levels and thus nonlinear lengths conveniently, but this is beyond the scope of the 
present work. 

The spectrum at the end of the propagation is shown in figure [7] We notice two main peak 
pairs appearing beyond Raman-gain bands (the two broad asymmetric peaks at Af = ±13 THz). 
The first pair is located at Af = ±34.8 THz, the second pair at Af = ±56 THz. The bandwidth 
of each individual peak belonging to the first pair is 0.25 THz, and the gain agrees well with our 
theoretical prediction (g w 100 km -1 ). The second PR peak pair exhibits smaller gain (g' w 70 
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Fig. 7. Output spectrum after the propagation in a periodically tapered PCF. All fiber pa- 
rameters are listed in Tab.[3] We detect the two Raman gain bands at A/ = ±13 THz, Stokes 
(red-detuned) and anti-Stokes (blue-detuned), the former exhibits as expected a larger gain. 
We also label the main two PR instability peak pairs which correspond respectively to the 
design requirement of Af = 35 THz and to the additional phase matching allowed for by 
FOD. 

km -1 ) and narrower bandwidth, and can be ascribed to FOD. which allows for an additional 
solution of the nonlinear phase matching condition, see e.g. l37l . In summary, the presence 
of j34 < has the advantage of leading to a period Tz larger than that predicted neglecting all 
HOD terms. It has also a drawback that an additional PR band appears, but we verified this 
applies, for our specific choice of parameters, only at the first PR order. If we include losses in 
our simulations we obtain a smaller available power, which affects mainly the peak gain, see 
Eq. (fTTt . Indeed the resonant detuning, see Eq. ( II Ob . is nearly independent on input power for 
a ~ 200, which corresponds to our a. 

5. Conclusion 

In this work we have presented a thorough analysis of parametric resonance instabilities occur- 
ring in a generalized nonlinear Schrodinger equations with varying dispersion and nonlinearity. 
We have shown that, in the case of GVD and nonlinear coefficients varying in a simple sinu- 
soidal way, it is possible to predict the maximum gain and instability margins. The calculation 
of the resonant frequencies is possible even in more general cases, since it depends only on 
the average values of coefficients and the period of variation. Our calculations provide analyti- 
cal estimates and are more accurate than those found in previous literature. We have validated 
our theory by means of numerical dynamical simulations and found a very good agreement. 
Finally we have designed a periodically tapered photonic crystal fiber which allows to achieve 
instability peaks at a large tunable frequency detuning from a given pump wavelength. Several 
higher-order resonant peaks are present but their gain and bandwidth are generally smaller than 
the one occurring at the smallest detuning. Such a system can be useful in quantum optical 
applications such as the efficient generation of entangled photon pairs in regions of frequencies 
far from the Raman gain peak. 
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